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Cosmological Simulations using GCMHD+ 
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ABSTRACT 

Radio observations of galaxy clusters show that the intra cluster medium is perme- 
ated by /iG magnetic fields. The origin and evolution of these cosmological magnetic 
fields is currently not well understood and so their impact on the dynamics of struc- 
ture formation is not known. Numerical simulations are required to gain a greater 
understanding and produce predictions for the next generation of radio telescopes. 
We present the galactic chemodynamics smoothed particle magnetohydrodynamic 
(SPMHD) code (GCMHD+), which is an MHD implementation for the cosmological 
smoothed particle hydrodynamic code GCD+. The results of 1, 2 and 3 dimensional 
tests are presented and the performance of the code is shown relative to the ATHENA 
grid code. GCMHD+ shows good agreement with the reference solutions produced by 
ATHENA. The code is then used to simulate the formation of a galaxy cluster with 
a simple primordial magnetic field embedded in the gas. A homogeneous seed field of 
10~^^G is amplified by a factor of 10^ during the formation of the cluster. The results 
show good agreement with the profiles found in other magnetic cluster simulations of 
similar resolution. 
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1 INTRODUCTION 



Cosmological magnetic fields are thought to be ubiqui- 
tous through out the universe and have been detected 
by radio observations on scales as large as galaxy clus- 
ters. On the galactic cluster scale, magnetic fields have 
been detected via Faraday rotation measurements of dis- 
tant quasars/ AGN, by the observation of radio galaxies 



chrotron emission from radio haloes (e.p 


\ Kronbere Il994l: 


Carilli k Tavlor 2002!: iGovoni k Ferettil 12004: IPizzo et all 


201ll: IClarke et all l200lh. The current 


observations find 



magnetic fields of /xG strength permeates the intra-cluster 
medium (ICM) of most galaxy clusters. Beyond the cluster 
scale, measurements of the magnetic field are far less cer- 
tain. Cosmological fields can be generat ed via Weibel' s insta- 
bility (Schlickeis er Shuklal 120031: iMedvedev et al.1 l2006ll 
Biermann's battery jBiermann^ *195rt ^Su bramanian et al.l 
ll994l : lKulsr"[^ et al E997; Gnedin et al. 20001)^ structure for- 
mation (^Harrison* ^70; Ichiki et al. 2006), galactic wind s 
(|Beck et al. 1996; Bertone et al. 2006; Donne rt et alll2009h . 
relativist ic charged particles (|Miniati k Bell 2011) and var- 
ious processes in the early universe fi vidrow. 2002 ). Cur- 
rent observations of the structural detail of these fields is 
limited, but the next generation of radio telescopes, such 
as the Square Kilom eter Array (SKA) (e.g. Gaensler 2006; 
I Johnston et al.ll2008l ) , will produce a wealth of observational 
data on the strength and structure of cosmological mag- 
netic fields. To compare these observations with our knowl- 
edge of cosmological magnetic fields we require numerical 



simulations. The predictions made by simulations allow for 
a comparison of the theory with the observation. This will 
produce a more detailed understanding of cosmological mag- 
netic fields. 

The evolution of primordial magnetic fields in 
a cosmological setting has been studied using both 
smoothed particle magnetohydr odynamics (S PMHD, e.g. 
Dolag Stasvszvn 2009; Bona fede et al. [2011>) and adap- 
tive mesh refinement c odes (e.g. iDubois Tevssierl l2008l : 
iMiniati Martini 1201 ll ). The difi"erent techniques produce 
good agreement on the predicted strength and profile of a 
magnetic field in a galaxy cluster. They show that the ini- 
tial structure of the field has little influence on the flnal fleld. 
The central cluster fleld strength predicted by these simu- 
lations agrees with the values inferred from observation. A 
few simulations have followed the creation of a cosmolog- 
ical magn etic fleld from eithe r galactic wind pollution of 
the ICM (iDonnert et al.l[2QQ9l ). AGN pollution of the clus- 
ter (|Xu et al.l I2OIOI ) or the Biermann batte ry mechanism 
during structure formation (IRvu et al .11 19 98^ to generate a 
fleld from zero fleld initial conditions. The magnetic fleld 
strength can also be predicted a posteriori using the veloc- 
ity flel ds from a purely hydrodynamical cosmological simula- 
tions (jRvu et al.ll2008l ). These simulations produce different 
fleld strengths than the values found from following the evo- 
lution of a primordial magnetic fleld, especially in fllament 
structures. Further simulations are required to examine the 
different origins of a magnetic fleld in the ICM and to test 
the validity of the predictions made. 

We present the implementation of magnetohydrody- 
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namics (MHD ) for the existing galactic chemo dynamics 
code (GCD+) (jKawata k Gibsonll2003l : iKawata et al.ll2009l . 
l201lh . GCD+ is a three-dimensional tree N-body/Smoothed 
Particle Hydrodynamics (SPH) code which incorporates self- 
gravity, hydrodynamics, radiative cooling, star formation, 
supernova feedback and metal enrichment. The following 
MHD implementation is fully compatible with all of the 
original features of the code and allows for a complete cos- 
mological simulation to be run. We discuss the choice of 
method to ensure tensile instability is avoided to an appro- 
priate level. We also discuss th e effect of an applied time 
dependant dissipation scheme ([Morris MonaghanI Il997l ) 
for the magnetic field and the effect of the Balsara switch 
(|Balsara ,1995) on the test simulations, including its need in 
a cosmological simulation. For the tests shown in this pa- 
per we ignore the additional processes and concentrate on 
the effect of introducing MHD to the simulations and the 
accuracy of the non-radiative solutions produced. As such, 
radiative cooling, star formation, supernovae feedback and 
metal enrichment are switched off in all results presented. 

In Section 2 we present the numerical implementation 
of the MHD equations used in GCMHD+. Section 3 shows 
the performance of the code in various test simulations. Sec- 
tion 4 shows the results of a cosmological simulation for the 
formation of a galaxy cluster with a simple primordial mag- 
netic field embedded in the gas. Our conclusions are then 
given in Section 5. 



2 MHD IMPLEMENTATION 

The hydrodynamic comp onents in the code are adopted 
from iKawata et al.l (|201lf ) . In this section we summarise the 
addition of the MHD implementation and the updated pa- 
rameters for the artificial viscosity for MHD simulations. We 
note that while our scheme does not ensure the V • B = 
constraint, we present results in Section 3.5 and 4 which 
show that it is well satisfied for all simulations. 



where Sf — 0.05/ii|V^iti|/y^ is the source term. Unless 
explicitly stated, A.C. is t urned on in aU simulations. 

The Balsara switch (Balsara 1995) reduces the A.V. 
when the code detects a shear flow. This prevents the ar- 
tificial viscosity from becoming the dominant force and gen- 
erating spurious forces in a shear flow. While running the 
MHD test suite with the code, especially shocktube test 5A 
in Section 3.1, it was found that this switch was causing 
the velocity to being captured poorly. We then removed the 
switch from the code and found that it produced a significant 
improvement in the code's ability to produce the velocity so- 
lution, with minimal negative effects on other tests. 

2.2 Induction Equation 

The magnetic field is evolved via the induction equation 
which, neglecting any form of dissipation and enforcing the 
V • B = constraint, takes the form 

IdB 

c dt 



(B- V)v-B(V-v) 



(2) 



This is the standard form of the induction equation and it is 
the correct choice as it is unaffected by magnetic monopoles. 
We then convert this to SPH components and the equation 
takes the form of: 



1 dB^ 
c dt 



1 



N 



t k Tjl Tjk I \ dWi Tij 



(3) 



where B are the components of the magnetic field and we 
sum over the / components, pi is the density at particle i, rrij 
is the mass of particle j, Vij is the velocity between the two 
particles, rij is the distance between the two particles and 
1/r^i is the factor correcting for the use of va riable particle 
smoothing lengths (jPrice Monaghanll2004bl l . The discrete 
form of equation ([2|) no longer enforces the V • B = con- 
straint. The magnetic field is allowed to act back on the fluid 
via a Lorentz force term in the momentum equation. 



2.1 Hydrodynamic Parameters 

Artificial viscosity (A.V.) is used in hydrodynamic codes 
to smooth discontinuities, such as shocks, in the velocity 
field of the simulation, thus allowing the code to capture 
the physics correctly. Applying a constant level of A.V. to 
the simulation cause s smo othing when it is not required. 
[m orris MonaghanI (|l997l l suggest an A.V. switch where 
the level of A.V. is allowed to vary in space and time for 
particles in the simulation. The minimum level of A.V. in a 
simulation is set by the parameter cind for a purely 

hydrodynamic simulat ion this is set to — 0.1(e.g, 



iRosswog Pried I2OO7I ). While running our MHD test suite 
for the MHD implementation we found some post shock os- 
cillation in the density field, hence for MHD simulations we 



increased the minimum A.V. to a„ 



0.6. 



We also i mplement artificial th ermal conductivity 
(A.C.) following iRosswog Pric3 (l2007l ). We allow the ther- 
mal conductivity parameter, , to vary with time for each 
particle. The parameter varies such that ^ ^ 2 via 



da? 
dt 



c 

Ti 



(1) 



2.3 Lorentz Force 



T he conservative form of the m agnetic stress tensor, derived 
bv lPhillips MonaghanI (|l985l \ is given by 



M 



47r V 



|2 ckl 
] 



)■ 



(4) 



The form of the stress tensor ensures conservation of mo- 
mentum at shocks. This generates a force that produces an 
additional term in the momentum equation, which in com- 
ponent form is 



dvl 
dt 



(mag) 



N 

47r ^ 



+ 



VLipl du \Yij\ 

Mf dWj r\j 
QjP^- du \rij\ 



(5) 



When the magnetic field dominates the gas pressure, this 
exactly momentum conserving form of the force becomes 
unstable. The magnetic stress can become negative and this 
leads to the tensile instability. The negative stress, or an 
attractive force, between particles can cause them to clump 
together and the simulation becomes highly unstable. This 
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requires a stable MHD code to have an additional term in the 
momentum equation to suppress the possibility of clumping, 
which takes the form of an instability correction. 

2.4 Tensile Instability Correction 

The Tensile Instability is a well known problem for SPMHD 
an d several m ethods have been proposed to suppress 
it. iMonaghanI (|2000) suggested the addition of an anti- 
dumping term to the momentum equation which prevents 
the occurrence of the tensile instability. At short distances 
a steepened kernel ensures that this additional term be- 
comes significant and the particles repel each other. This 
was found to be effective for one and two dimensional simu- 
lations. However, for 3D simulations with variable smooth- 
ing lengths this method is no longer effective at preventing 
clu mping. 



lB0rve et al.l (|200lh suggested directly subtracting any 
non-vanishing monopole terms from the momentum equa- 
tion. This can be implemented by an additional term in the 
momentum equation, which takes the form 



dt 



k \ (corr) 



-—Bi > rrii 
47r 



_B^dW^r\ 



+ 



B 



du 



(6) 



where the parameter /3 controls the level of non-vanishing 
divergence subtracted. In principle the the addition of this 
term to the momentum equation breaks the momentum con- 
servation of the formulation, but it does not seem to cause 
any major effects in our te st simulations. lB0rve et al.l (| 200 il l 
used a value of ^ = 1 and IPolag Stasvszvnl (|2009l ) found 
from testing that p = 1 produced no harm to their results. 
lB0rve et all ('2004') suggested that stability can be achieved 
with /3 < 1 and that /S = 0.5 should be used to minimise 
any non- conservative contribution. After testing we find that 
13 = 0.5 produces the optimal results for our code and that 
it removes the tensile instability effectively (see Section 3.5). 



2.5 Artificial Magnetic Dissipation 

When no magnetic field is present a good estimate of the 
speed at which a signal can be sent from one particle to 
another is given by 

= ^(^^ + ^■?') ~ f^^^^j ■ ^'^j- c^) 

where v^j^ is the signal velocity between the particles i and 
j, a is the sound speed at particle i and Vij is the velocity 
between the particles. In the presence of a magnetic field a 
variety of MHD waves can propagate. The simplest gener- 
alization of the signal velocity is to replace the sound speed 
with the fastest magnetic wave (Pr ice Monaghan. .200431 ). 
The sound speed is then given by 



1 

7^ 



c? + 



47rpi 



(8) 



In order to treat MHD shock fronts correctly a dissipation 
term is required to resolve any steep gradients in the mag- 
netic field. We implement an artificial magnetic dissipation 
analogous to an artificial viscosity, ba sed on the change 
of the total magnetic field, following IPrice MonaghanI 
(|2004al l. This artificial dissipation is included via an addi- 
tional term in the induction equation, which takes the form 



c\ dt J 2 

x(B, -B,)^ • V^m. 



(9) 



The strength of the dissipation is controlled by the parame- 
ter . The dissipation will lead to a generation of entropy 
according to 



(dai\ 
\~dtJ 



7 - 1 

N siq 

i=i 



2 ^i: 



V^T4^^,,(10) 



where a = (P/p^) (jSpringel Hernauistll2002h . It was found 
that the artificial dissipation significantly reduced the noise, 
however with a constant val ue of it also lead t o a sm ooth- 
ing out of acute features. iPrice Monagha 3 (|2005l l pro- 
posed that this could be avoided by making indepen- 
dent and time varying for each particle. We allowed to 
vary by integrating an equation very similar to the one for 
time-dependant viscosity. 



da 



[aB — Q^ri 



(11) 




where a^*^ is the minimum level of dissipation applied to 
each particle in the simulation and the source term, 5*, was 
chosen to be 



(12) 



The parameter r controls how fast the dissipation decays. 
This combined with the signal velocity defines the distance 
from the shock for the artificial dissipation to return to the 
minimum value. It is defined by 

hi 



where hi is the smoothing length for particle i and a value 
of 0.2 was chosen for the constant C, i.e. 5 smoothing 
lengths. In order to conserve momentum the average value 
— \{af + af) is used in all simulations. Tests indi- 
cate that the use of a particle variable dissipation parame- 
ter greatly improves the code's shock capturing capabilities 
without significant smoothing of sharp features. 



3 TEST SIMULATIONS 

Having outlined the additions made to produce the 
GCMHD+ code, we run a series of test simulations which 
have become the standard set to show the performance of a 
numerical MHD scheme. We ran a large number of series of 
test simulations with different parameters associated to the 
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Test 


Nl 


PL 


Vl 


Bl 


Pl 


Nr 


PR 


Vr 


Br 


Pr 


lA 


540 


1.00 


(10.0,0.0,0.0) 


(5.0,5.0,0.0)/(47r)0-5 


20.0 


540 


1.00 


(-10.0,0.0,0.0) 


(5.0,5.0,0.0)/(47r)0-5 


1.00 


IB 


1000 


1.00 


(0.0,0.0,0.0) 


(3.0,5.0,0.0)/(47r)0-5 


1.00 


100 


0.10 


(0.0,0.0,0.0) 


(3.0,2.0,0.0)/(47r)0-5 


10.0 


2A 


540 


1.08 


(1.2,0.01,0.5) 


(2.0,3.6,2.0)/(47r)0-5 


0.95 


500 


1.0 


(0.0,0.0,0.0) 


(2.0,4.0,2.0)/(47r)0-5 


1.00 


2B 


1000 


1.00 


(0.0,0.0,0.0) 


(3.0,6.0,0.0)/(47r)0-5 


1.00 


100 


0.10 


(0.0,2.0,1.0) 


(3.0,1. 0,0.0)/(47r)0-^ 


10.0 


3A 


550 


0.10 


(50.0,0.0,0.0) 


-(0.0,1. 0,2.0)/(47r)0-5 


0.40 


550 


0.10 


(0.0,0.0,0.0) 


(0.0,1.0,2.0)/(47r)0-^ 


0.20 


3B 


550 


1.00 


(-1.0,0.0,0.0) 


(0.0,1.0,0.0) 


1.00 


550 


1.00 


(1.0,0.0,0.0) 


(0.0,1.0,0.0) 


1.00 


4A 


1000 


1.00 


(0.0,0.0,0.0) 


(1.0,1.0,0.0) 


1.00 


200 


0.20 


(0.0,0.0,0.0) 


(1.0,0.0,0.0) 


0.10 


4B 


400 


0.40 


(-0.6699,0.9826,0.0) 


(1.3,0.0025,0.0) 


0.5247 


1000 


1.000 


(0.0,0.0,0.0) 


(1.3,1.0,0.0) 


1.00 


4C 


650 


0.65 


(0.667,-0.257,0.0) 


(0.75,0.55,0.0) 


0.50 


1000 


1.000 


(0.4,-0.94,0.0) 


(0.75,0.0,0.0) 


0.75 


4D 


1000 


1.00 


(0.0,0.0,0.0) 


(0.7,0.0,0.0) 


1.00 


300 


0.300 


(0.0,0.0,1.0) 


(0.7,1.0,0.0) 


0.20 


5A 


960 


1.00 


(0.0,0.0,0.0) 


(0.75,1.0,0.0) 


1.00 


120 


0.125 


(0.0,0.0,0.0) 


(0.75,-1.0,0.0) 


0.10 



Table 1. Summary of the initial conditions for all ID test simulations, where is the number of particles, p is the density, V is the 3D 
velocity structure, B is the 3D magnetic field structure and P is the pressure. L and R denote the left and right halves of the simulation. 
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0.2 
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Figure 1. Result for the shocktube test lA. The ATHENA reference solution is shown by the red line and the blue points show the 
result produce by the code. They are in general in good agreement except for a small amount of deviation in the density and magnetic 
profiles. 



MHD implementation, such as a^i^ ^^id /3. This lead to our 
best parameter set for the suite of test simulations. This set 
is: a^ir, 0.05, a^ax 1-0, ^ 0.5, a^ri 0.6, artificial 
conductivity on and the Balsara switch turned off (see Sec- 
tion 3.5). However, as demonstrated in Section 4 for the cos- 
mological simulation we require a^^^i = 0.0. Therefore, we 
set a'^in — 0.0 as a fiducial case. In this section we present 
the test simulation suite results with this fiducial set of pa- 
rameters ( a^i^ 0.0, a^^^ 1.0, P 0.5, 0.6, 
artificial conductivity on and the Balsara switch turned off) 
and demonstrate that o;^^^ = 0.0 leads to a satisfactory 
result for all the test simulations. Due to the complex na- 
ture of MHD interactions no analytical solution exists for 
many of the tests presented below and so the performance 
of the code is compared to a reference solution. This is pro- 



yided by the publ icly available ATHENA MHD mesh code 
(jStone et al.ll2008l ). 



3.1 Magnetic Shock Tubes 

Magnetic shocktubes are the standard test of any numer- 
ical MHD scheme. The addition of MHD allows for more 
complex solutions to the Riemann problem due to the pres- 
ence of slow and fast Alfven waves. Sharp magnetic features 
present in the simulation will be smoothed by the artificial 
dissipation and the full effects of any applied regularization 
scheme must be explored. This means that a f ull set of shock 
tube problems, as presented in IRvu JonesI (1995), are re- 
quired to rigorously test the effects of the applied artificial 
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-0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 

X X 



Figure 2. Result for test IB. The code shows good agreement with the reference solution in the density and pressure profiles. The 
velocity profile is not captured exactly in the low density region. There is a small amount of noise present in the magnetic field solution. 




-0.4 -0.2 0^0 0.2 0.4 -0.4 -0.2 0^0 0.2 0.4 



Figure 3. Result for shocktube test 2A. The two solutions show good general agreement except for a small amount of noise present in 
the density, pressure and magnetic field. 



magnetic dissipation scheme and to establish the best pa- 
rameters for ID solutions produced by the code. 

The set up of all the one dimensional shocktube sim- 
ulations can be found in Table [1] While the particle posi- 
tion is only allowed to vary in the x direction, the parti- 
cle velocity and the magnetic field are allowed to vary in 



3 dimensions. As no analytic solution is known, the code 
is compared against a reference solution obtained using the 
ATHENA code. The resolution of the ATHENA simulation 
depended on the shock tube being run, but the parameters 
of each test were left unchanged from the originals provided. 
Figs. 1-11 show the density, pressure, Vx and By out- 
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Figure 4. Result for shocktube test 2B. There is a small amount of noise in the magnetic field profile and the code doesn't quite capture 
the velocity exactly. 
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Figure 5. Result for the shocktube test 3A. The pressure, velocity and magnetic field are well captured by the code when compared to 
the reference solution. Due to the difficult hydrodynamic set up the density is not well captured for the transition. 



puts for all the shocktube tests. The code captures the ma- 
jority of the features precisely. For the tests, lA and IB, the 
code captures all three stages, including the intermediate 
phase, for both the strong and weak shock accurately. For 
the density, at x = 0.02, and magnetic field, at x = 0.00 
and X = 0.02, of lA and the velocity, at x = 0.39, and mag- 



netic field, at X = —0.14, of IB there are small deviations 
from the reference solutions. The code's ability to capture 
the velocity in test IB improves as the resolution increases. 

The second set of tests, 2 A and 2B, have a 3D velocity 
and magnetic field structure. All the features are well cap- 
tured by the code and only the transition in the magnetic 
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-0.4 -0.2 0^0 0.2 0.4 -0.4 -0.2 0^0 0.2 0.4 



Figure 6. Result for shocktube test 3B. The code captures all of the variables with very little noise. There is a discrepancy between the 
code and the reference at x = 0.0 for the magnetic field. 




-0.4 -0.2 0^0 0.2 0.4 -0.4 -0.2 0^0 0.2 0.4 



Figure 7. Result for shocktube test 4A. The density, pressure, velocity and magnetic field profiles all agree with the solution from 
ATHENA. There is a small amount of noise in the velocity and magnetic field profile. 



field of test 2B, between the intermediate and lower state, 
at X — 0.05, shows any visible deviation from the reference. 
In both tests there is also a small amount of oscillation vis- 
ible in the magnetic field at the shock boundaries. Test 2 A 
shows some oscillation in the density and pressure as well. 
The third tests, 3 A and 3B, show the code's ability to 



handle magneto-sonic features. The pressure, velocity and 
magnetic field profiles produced by the code for test 3A agree 
well with the reference solution. The density profile in test 
3 A is poorly captured. There is a large dip in the density 
compared to the reference solution. This is due to the hard 
hydrodynamical conditions for the test, where some of the 
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-0.4 -0.2 0^0 0.2 0.4 -0.4 -0.2 0^0 0.2 0.4 



Figure 8. Result for shocktube test 4B. The result produced by the code shows good agreement with the reference solution for all of 
the parameters. The is some noise where the magnetic field switches off. 
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Figure 9. Result for shocktube test 4C. The density, pressure and velocity agree with the reference solution. The magnetic field shows 
good agreement except at the shock where there is significant noise present. This can be reduced by using the test optimized parameters. 



particles are colliding at high speed with stationary ones, 
producing a wall heating error. This is effectively a wall for 
the high speed particles and due to the restriction of only 
being allowed to travel in one spatial dimension they crash 
in to this wall and oscillate back and forth. This causes the 
resulting dip in the density at the wall. This problem is not 



experienced by the mesh code and so it produces the cor- 
rect solution. The velocity profile for test 3B agrees with 
the solution produced by ATHENA. For the density, pres- 
sure and magnetic field solutions are well captured with a 
small deviation at x = 0.0. This code shows smaller peaks 
compared to the ATHENA result. However, the result from 
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-0.4 -0.2 0^0 0.2 0.4 -0.4 -0.2 0^0 0.2 0.4 



Figure 10. Result for shocktube 4D. The density, pressure, velocity and magnetic field all show good agreement with the reference 
solution. There is a small amount of noise present in the magnetic field solution. 




-0.4 -0.2 0^0 0.2 0.4 -0.4 -0.2 0^0 0.2 0.4 

Figure 11. Result for shocktube 5A. The code produces a solution which agrees with the reference for the density, pressure and magnetic 
field with very little noise present. The velocity is poorly captured in the low density region. 



IRvu Jone j (l995") for test 3B is more comparable to the 
solution produced by GCMHD+. We therefore believe the 
inconsistencies between the results could be due to an arti- 
fact in the ATHENA solution. 

The fourth set of tests deal with features caused by 
the magnetic field turning on and off behind fast and slow 



moving rarefactions and shocks. The code produces density, 
pressure and velocity profiles which show good agreement 
with the reference solutions. The magnetic field is well cap- 
tured, with a small amount of oscillation, in 4A, 4B and 4D. 
The magnetic field profile found in test 4C does agree with 
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using ATHENA (red line). The code captures the majority of the 
density features very well. There is a small difference between the 
ATHENA solution and the code either side of the main density 
peaks due to smoothing. 




Figure 12. Plots of the density (upper) and the magnetic pres- 
sure (lower) for the fast rotor test at t = 0.1. The arrows show the 
size and direction of the magnetic field. The shape and position of 
the main features agree with other rotor test solutions published 
in the literature. 



the reference solution, but there is a very large amount of 
oscillation present in this test at the shock front. 

The final test, 5 A, is the commonly used shock tube 
test of iBrio Wul (p,988). This test involves a shock and 
rarefaction moving together. The solution produced by the 
code captures all of the features very well for the density, 
pressure and magnetic field, with a small amount of oscilla- 
tion visible in all of the profiles. The velocity is well captured 
in the high density region, but it is poorly captured in the 
low density region. Overall, the code captures all of the fea- 
tures present in the tests and produces very similar results 
to the reference solutions provided by ATHENA. 



3.2 Rotor test 

The fast rotor test has been used many times to check th e 
validity of solutions produced an MHD code (jTothI 
The 2D test is set up with a dense rotating disc embedded in 
a low density, static ambient medium. A uniform magnetic 
field is applied across the entire simulation. A constant field 
is set in the x direction with a strength of Bx — 2.5/v^- The 
rotating disc has a radius ro = 0.1, initial density of p = 10 
and pressure of P = 1. The rotational velocity is given by 
Vx = 2{y - 0.5)/ro, Vy = 2(x - 0.5)/ro and Vz = 0. The 
ambient medium has a density p = 1 and pressure P = 1. 

One way to produce the density contrast is to use par- 
ticles of different mass for the disc and the background 
medium. However, this produces spurious unwanted effects. 
Instead we apply the same mass to all the particles and put 
ten times more particles in the disc region. First, the am- 
bient medium is laid down in a regular hexagonal lattice. 
Then a second lattice is placed with smaller particle spacing 
in a region 2ro x 2ro centred at (0, 0). The larger separation 
particles for which r < ro and smaller separation particles 
with r > ro are removed. This ensures that all particles in 
the simulation have the same mass. A hexagonal lattice is 
used instead a standard square lattice to reduce any discon- 
tinuities between the disc and the background. The lattice 
results in a set up of 400 x 460 particles in the background. 
This results in a total of 236626 particles in the simulation. 

The simulation is evolved until t = 0.1 and the result 
can be qualitatively seen in Fig. 1121 The material contained 
in the disc is thrown out into the surrounding medium, but 
is contained by the magnetic field. The quantitative com- 
parison between the code and the ATHENA simulation can 
be seen in Fig. 1131 The ATHENA simulation was run with 
400 X 400 cells. The result agrees well with the solution pro- 
vided by ATHENA and shows very little smoothing of fea- 
tures by the time dependant artificial dissipation. There is a 
small deviation between the two solutions either side of the 
main density peaks. The difference between the solutions is 
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Figure 14. The density (upper row) and magnetic pressure (lower row) distributions for the Orszag-Tang vortex at t = 0.5 for three 
resolutions: 128 x 146 (left), 256 X 294 (centre) and 512 x 590 (right). The magnetic field direction and strength is shown by the arrows. 
The initial velocity field and periodic boundaries lead to complex interactions between the magnetic field and shock present in the 
simulation. The code produces results very similar to others presented in the literature. 



caused by the code smoothing out the density change across 
the edge of the shock front. 



3.3 Orszag-Tang Vortex 

The compressible Orszag;-Tang: vortex w as developed from 
a test problem in lOrszag Tand (| 19791 ) and is a common 
test of MHD implementations. It shows the code's ability 
to handle the interaction between different classes of shock 
waves and the transition to MHD turbulence. Using 7 = 5/3, 
a magnetic to thermal pressure ratio of 10/3 and an aver- 
age Mach number of unity, a uniform medium with peri- 
odic boundaries is set out in a hexagonal lattice with P = 
5/3^0 aiid p — 7P/V0, where Bq — 1/V^ and vq — 1.0. 
The magnetic field is initially set as Bx — — Bo sin(27r^). 
By = Bq sin(47rx) and Bz = 0. The gas is also given an 
initial velocity of Vx = — sin(27r^), Vy = i;o sin(27rx) and 
Vz — 0. We performed the simulation for 3 different resolu- 
tions: 128 X 146, 256 x 294 and 512 x 590. 

The simulation was evolved until t = 0.5 and the re- 
sults are shown in Fig. 1141 The results show go od agreement 
with t hose published in the literatures, such as IRvu JonesI 
(|l995l ). As the resolution increases the complex interplay be- 



tween the four shock fronts and the magnetic field becomes 
clearer, with more small scale features easily visible in the 
512 X 590 run compared to the 128 x 146 run. 

For a quantitative comparison of the code's ability, we 
measure the variation of the pressure as a function of x at 
a fixed y and compare the results to the solution produced 
by ATHENA. The results are seen in Fig. [15] There is good 
agreement between the two solutions. A small amount of 
smoothing can be seen in the upper plot at x = —0.45 where 
the ATHENA solution produces very sharp features. 



3.4 MHD Point-like Explosion 

The MHD blast wave test is identical to the Sedov test com- 
monly used for testing pure-hydrodynamic codes but with a 
magnetic field added to the simulation. The uniform field is 
initially added in the x direction, such that B = [Bo, 0,0], 
where Bo = 3. The medium is a constant density, p = 1, 
with a hot point source embedded in it. The simulation box 
runs from —0.5 to 0.5 in the x, y and z directions. The 100^ 
particles are set out on a cubic grid ensuring that a parti- 
cle occupies the position at (0,0,0). The central particle is 
then given an energy 100 times the energy of the ambient 
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Figure 15. The pressure along ID cuts dX y = 0.3125 (upper) 
and y = 0.4277 (lower) in the Orszag-Tang vortex. These two cuts 
have been chosen to allow comparison with other Orszag-Tang 
solutions presented in the literature. The ATHENA solution is 
plotted in red, while the solution produced by the code is shown 
via the blue dots. There is good agreement between the solutions. 

medium. We do not apply the smooth ed central hi^h en- 
ergy sphere used in some literature (e.g. IPolag Stasvszvnl 
I2OO9I ). Instead we let the energy spread from a single high 
energy particle. This is a much harder initial condition than 
the smoothed case. Due to the strength of the field, the mag- 
netic pressure plays an important role in the evolution of the 
shock. 

As Fig. [16] shows, particles moving in a direction per- 
pendicular to the orientation of the field are constrained by 
the magnetic tension force, which prevents the density at 
the shock front from increasing to the levels seen in the x 
direction. It is also noticeable that the field lines are not sig- 
nificantly bent by the shock front. The solution produced by 
the code agrees well with ot her 3D MHD blast wave results 
shown in the literature fe.g. lBalsarall200ll : iRosswog &: Pried 
I2OO7I ). 

3.5 Numerical Parameters 

The simulations shown in previous sections use our best 
compromised set of numerical parameters that can produce 
satisfactory results for all of the tests, including the cosmo- 
logical simulation. The value of the magnetic dissipation is 




Figure 16. Plot of the density (top) and pressure (bottom) for 
the 3D magnetic blast wave test. The arrows indicate the strength 
and direction of the magnetic field in the x-y plane. The solution 
produced by the code agrees well with other solutions published 



a compromise between minimising the smoothing of sharp 
features and reducing oscillations in the magnetic field. This 
can be seen in the cosmological simulation, where the mag- 
netic field is highly dependant on the applied dissipation, 
and other simulations which show the presence of oscilla- 
tions in the solution, such as the magnetic field for test 4C 
in Section 3.1. In order to get more of a quantitative mea- 
sure of the quality of the solution the mean of all V • B / 
errors, Ss, in the simulation volume can be calculated using 
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Figure 17. Mean divergence error as a function of the minimum 
applied magnetic dissipation, o;^^^. The error is averaged across 
all particles in the simulation. The labels lA to 5A represent the 
ID shock tube tests in Section 3.1, ROT is the 2D rotor test 
in Section 3.2, OTl and OT2 are the Orszag-Tang vortex at a 
resolution of 128 x 146 and 256 x 294 respectively (Section 3.3) 
and MBW is the MHD Point-like explosion test (Section3.4). 

The numerical parameters can then be varied and the er- 
ror can be plotted as a function of the parameter. Fig. [TTl 
shows the mean error of the simulation as a function of the 
minimum applied artificial dissipation, a^^^ for all of the 
test simulations shown in Sections 3.1-3.4. Fig. [iTl shows 
that applying a small amount of dissipation to the simu- 
lation improves the divergence error, or the noise, of the 
solution. However, it also shows that increasing a-^i^ be- 
yond Oi'^^n = 0.1 does not lead to a further reduction in the 
noise, even in multi-dimensional tests. Any small amount of 
applied dissipation will lead to some smoothing of features 
and so the best range for for the non-cosmological test 
simulations was to allow it to vary between 0.05 and 1.0 for 
each particle. However, due to the dependence of the cosmo- 
logical magnetic field on the level of applied dissipation the 
value of was allowed to vary between 0.0 and 1.0, and 
therefore we chose a-^i^ — O-O- Fig [T3 demonstrates that 
this choice of a'^^^ provides a reasonably low level of error. 
The second important parameter for GCMHD+ is 
which controls th e level of V • B force subtraction. 

In lB0rve et al. I iooi) 

it was suggested that the value of 
P could be less than on e and stability still ensured. In 
iDolag k, Stasvszvnl (|2009l l they argue that it is unclear as 
to whether the conclusions by B0rve et al. hold in 3D. We 
ran all of our test simulations and measured the value of 
as a function of ^. Fig. [THl shows that increasing p has a dif- 
ferent effect depending on the dimensions of the simulation. 
Note that if we apply f5 < 0.2 the code collapses and can- 
not solve the test problems. For ID tests, an increase in p 
generally leads to an increase in the average error. However, 
for higher dimension tests an increase leads to a reduction 
in the error, but with a reduced return for a value of 0.5 
or more. From the result shown in Fig. [18] and to minimise 
the violation of momentum conservation, we chose a value 
of = 0.5. 




/3 

Figure 18. Mean divergence error as a function of /3. A value of 
/3 < 0.2 does not fully remove the tensile instability and the code 
does not properly run. The labels have the same meaning as Fig. 
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Table 2. Properties of the cluster simulations. All other param- 
eters were kept constant to allow the effect of dissipation and 
thermal conductivity to be seen. 

4 COSMOLOGICAL SIMULATION 

The Santa Barbara Cluster simulation was used as t h e base 
for the cosmological simulation. Since 'Frenk et all (|l999h 
this has become a standard test for cosmological hydro- 
dynamic codes. The initial conditions ar e set u p to create 
a galaxy cluster. It was shown in Fren k et al.1 |l999l) that 
different codes and numerical techniques produce generally 
consistent results. 

The simulation assumes a classical, flat CDM cosmol- 
ogy. The simulation allows a spherical region of 32/i~^ to 
expand with the Hubble flow. The test uses open boundary 
conditions. The gas particles have a mass of 8.67 x lO^M© 
and the dark matter particles have a mass of 7.80 X 10^ Mo . 
(This corresponds t o a slightly low er resolution than the 1 x 
simulation of Dola g Stasvszv"nl (j2009|)). To this a homo- 
geneous magnetic field of 10~^^G is initially applied. The 
simulation is started at a redshift of z = 20 and evolved 
to the current epoch. This produces a final cluster of mass 
1.16 X 10^^ M0. The cluster simulation was run varying the 
value of a-^in to test the effects of the dissipation scheme on 
the development of the magnetic field. The effect of artificial 
thermal conductivity (A.C.) on the magnetic field was also 
tested by running the simulation with and without it. The 
values for the different runs can be seen in Table [2] 

Fig. [19] shows the hydrodynamic properties of the sim- 
ulated clusters. The radial profiles of the density, tempera- 
ture and velocity dispersion for the simulated clusters with 
artificial conductivity agree well with each other. The clus- 
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Figure 19. Radial profiles of the density (top), temperature 
(middle) and velocity dispersion (bottom) for all of the simu- 
lated clusters at z = 0.0. These agree well with each other and 
with the results in the literature. 



ters without artificial conductivity show a higher core den- 
sity, lower core temperature compared with the c lusters that 
have artificial conductivity (|Kawata et al.|[201ll V but their 
profiles agree between models with different values of a^^^. 
This result is expected as the magnetic energy density is 
very much less than the kinetic energy and so any mag- 
netic field present should not influence the development of 
the cluster. The proflles of the clusters which use artifl- 
cial conductivity agree well with the mesh code results pre- 
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Figure 20. Radial profiles of the magnetic field, mean dissipation 
(|q;^|) and mean V • B error for all of the simulated clusters at 
z = 0.0. The magnetic field shows that enforcing a minimum 
limit threshold of dissipation causes the field to decay too quickly. 
The middle plot shows that the simulation will settle at a lower 
value of dissipation, and a higher level of magnetic field, when 
^min ~ applied. The bottom plot shows that V • B error is 

kept at a satisfactory level for all of the simulated clusters. 

sented in lPrenk et al with a core density peaking at 

10 ^^ M(:) pc~^ and teni perature of 10^ which is discussed 
in iKawata et alJ (|201lh . 

The magnetic field of the simulated cluster can be seen 
in the upper plot of Fig. 1201 Clearly the final profile and 
strength of the magnetic field depends on the minimum 
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value of dissipation enforced on the simulation. Clusters with 
higher levels of minimum dissipation have a lower strength 
magnetic field through the entire volume, and the strength 
of the field in the core is reduced by 4 orders of magnitude 
compared to the case with a^^^ — 0-0 (simulation CI). This 
strong dependence on the level of dissipat ion was also seen 
in results presented in IPolag Stasvszvnl (^009). The clus- 
ters simulated without artificial conductivity show a similar 
radial profile for their magnetic field when compared with 
the cluster simulated with A.C. and the same minimum dis- 
sipation, but the strength of the field is greater. We found 
that A.C. lead to simulations with a much smoother den - 
sity and temperature distribution fsee lKawata et al]l201ll ). 
Simulations without A.C. have more small scale structure 
in their density field, which amplifies the magnetic field and 
keeps it at a higher level compared to the models with A.C. 

The middle plot of Fig. [20] shows the radial profile of 
i.e. the average level of dissipation applied. It shows 
that enforcing even a very small level of minimum dissipa- 
tion is too high for a cluster simulation. The radial profile 
for CI shows that when a^^^ — 0-0 is applied, the code 
produces a dissipation of roughly 10~^. The average value 
IS a = 3.35 X lO"'^. The cluster core has a lower level 
of dissipation compared to the edges of the cluster, where 
a small amount of material is still falling in. The effect of 
these minor mergers can be seen in the top panel of Fig. [20] 
at roughly 1 Mpc. 

Recentlv lBonafede et al.l (| 20 111 ) used a constant dissipa- 
tion level of ?7m = 6 X lO^^cm^s"^ to produce the observed 
magnetic field in their cosmological simulation. They also 
derived a value of rfm — 2 x lO^^cm^s"^ from turbulence 
arguments. Their rim simply corresponds to /2. We find 
an average value of — 3.35 x 10~^. Converting this to cgs 
units, we obtain rfm — 2.14 x lO^^cm^s"^ which is similar to 
their values. 

The lower plot of Fig. [20] shows the radial profile for 
the average V • B error calculated by eq. (14). The error 
remains at the level of a few percent throughout the cluster 
volume and is acceptable for all levels of applied magnetic 
dissipation. 

The evolution of the magnetic field from z — 1.5 to 
the current epoch in simulation CI is shown in Fig. 1211 At 
z = 1.5 there are two major protoclusters and these then 
merge before z — 1.2. The peak field in the cluster remains 
roughly constant from z — 1.5 to the final value found in 
the simulation. The effect of merging and in-falling material 
can be seen as the radial field profile evolves and the change 
in field profile at z = 1.2 is due to the major merger of 
the two large protoclusters. The field then relaxes back via 
dissipation to its final shape. 

There is a very slight increase in the peak field at z — 
0.6. This is when the last big in-fall of material occurs for 
the cluster. This shows that in the absence any additional 
source terms the magnetic field is amplified by the merging 
of protoclusters and by minor mergers at later epochs. 

Fig.l22lshows the evolution of the overall magnetic field 
strength over the redshift range of z = 1.5 to z = 0.0. The 
overall magnetic field strength grows as the cluster assembles 
itself. This is due to the in-fall of material on to the cluster. 
This accretion causes amplification of the magnetic field, 
increasing its strength. The peak changes with redshift due 
to the increase in size of the cluster. As the cluster grows. 
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Figure 21. Radial profiles for the magnetic field for the redshift 
range of z = 1.5 to z = 0.0 of simulation CI. The magnetic field 
in the cluster has already roughly reached it maximum value and 
correct profile by the time the final cluster is forming at z=1.5. 
The in-fall of small structures causes the value to change slightly. 
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Figure 22. Radial profiles for the strength of the magnetic field 
for the redshift range of z = 1.5 to z = 0.0 of simulation CI. As 
the cluster evolves the strength of the magnetic field increases. 
The power peaks at increasing radius as the cluster grows. 

the largest possible coherence length of the magnetic field 
will increase accordingly and so the peak of the strength will 
tend to larger radii. In this simulation the strength profile 
at z = 1.2 is very different to the others due to the effects 
of the major merger disturbing it. 

The Fourier transform of the magnetic field was calcu- 
lated over the same redshift range as Fig. 21. The edge of 
the cluster was defined by its virial radius and this set the 
box size for each transform. We adopt the comoving scale, al- 
lowing the different redshift profiles to be directly compared. 
The Fourier power spectra of the magnetic fields are shown 
in Fig. 1231 The general trend shows an increase in the power 
of the magnetic field as the cluster evolves. At high frequen- 
cies the power increases very little as the cluster grows, but 
at lower frequencies the field increases between z = 1.5 and 
^ = 0.0 implying that the field becomes more coherent as 
the cluster evolves. The major merger is clearly visible as the 
low frequency power spectrum increases very rapidly as the 
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Figure 23. The power of the Fourier transform of the magnetic 
field as a function of frequency for simulation CI, between the 
redshift of z = 1.5 to z = 0.0. The plot shows a clear increase in 
the low frequency power of the field as the cluster evolves, while 
the higher frequency power increases very little. The dotted line 
represents a cut-off where window effects dominate the spectrum. 

two protoclusters merge between z = 1.5 and z — 1.2. (Note 
that any frequency lower than the dotted line is dominated 
by the window function and should be ignored.) 



5 SUMMARY 

We have introduced the MHD component to the N- 
body/SPH code GCD+. We discussed the addition of the 
equations of ideal MHD, the choice of instability correction 
and the addition of dissipative terms to treat discontinuities 
in the magnetic field. We implement schemes to re move the 
tensile instability, suggested bv fB^rve et al l (l200lh. and for 
artificia l magnetic dissipation, following IPrice &: Monaghanl 
The code's ability to vary this dissipation in the 
simulation and allow each particle to evolve it own dissi- 
pation constant is presented. We put the code through a 
set of standard ID shocktube tests, the fast rotor test, the 
Orszag-Tang vortex test and the MHD Point like explosion 
test. The numerical parameters were varied for all tests and 
the best compromise between noise reduction and minimised 
smoothing was found. The code with the best compromised 
parameter set performs very well in these tests and agrees 
with the reference solutions provided by the ATHENA mesh 
code, where they are available. The code shows no sign of 
the tensile instability and the magnetic dissipation scheme 
produces very little smoothing, while allowing the code to 
accurately capture the features. We then applied the code 
to a cosmological simulation for the formation of the Santa 
Barbara galaxy cluster. The code produces the expected hy- 
drodynamic parameters of the cluster. The magnetohydro- 
dynamic parameters are also well captured and show a sim- 
ilar leve l of magnetic field to the l iteratures with similar res- 
olution (Dolag & Stasyszyn 2009). We demonstrate that no 
minimum limit for the parameter of the dissipation of mag- 
netic field, a^in — O-Oi is necessary to minimise the artificial 
dissipation for the cluster scale magnetic field. This require- 
ment is significantly lo wer than the previous SPMHD imple- 
mentations, except for lBonafede et aP (|201lh . Our extensive 



test simulations in Section 3 demonstrate that cx^in — 0.0 
still leads to satisfactory results. Encouraged by the suc- 
cess of our new MHD code, GCMHD+, we will apply it to 
higher resolution cosmological simulations, and study how 
magnetic fields developed in the evolving universe. 
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APPENDIX A: CODE PARAMETERS 

Sections 3 and 4 show test results for the best compromised 
parameter set of GGMHD+. Ghanging these parameters will 
change the solutions produced by the code for each test. Here 
we present the effect of changing these parameters. 

In GGMHD+ the Balsara switch has been removed from 
the implementation. When using the switch we found that 
the velocity solutions produced for shocktube 5A were incor- 
rect. Fig. lAll shows the Vx solution produced by GGMHD+ 
with and without the Balsara switch. When the switch is in- 
cluded the velocity in negative velocity region is incorrectly 
captured. When the switch is removed the velocity solution 
produced by the code shows improved agreement with the 
reference solution provided by ATHENA. For this reason the 
switch was removed. 
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Figure Al. Plot of Vx for shocktube 5A. The reference solution, 
provided by ATHENA, is shown in red. Two solutions produced 
by GCMHD+ with (yellow diamonds) and without (blue crosses) 
the Balsara switch are shown for Vx ■ The solution in the negative 
velocity region is improved by removing the Balsara switch. 
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Figure A2. Plot of the density for shocktube 3A. The ATHENA 
solution is shown by the red solid line, the code's solution without 
A.C. is shown by the blue points. This is significantly worse when 
compared to Fig. 5. 
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Figure A3. Plot of Vy for shocktube 5 A. The reference solution, 
provided by ATHENA, is shown in red. The solution produced 
by GCMHD+ with /3 = 0.2 is shown by the blue circles. 
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Figure A4. Plot of Vy for shocktube 5A. The reference solution, 
provided by ATHENA, is shown in red. The solution produced 
by GCMHD+ with /3 = 0.5 is shown by the blue circles. 
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Figure A5. Plot of Vy for shocktube 5A. The reference solution, 
provided by ATHENA, is shown in red. The solution produced 
by GCMHD+ with /3 = 1.0 is shown by the blue circles. 

In Section 3.1 we showed the results for the ID shock- 
tube test 3A. The density profile produced by the code shows 
evidence of a wall heating error at x — 0.25. A wall heating 
error can be reduced by applying A.C. at the point where 
it occur s. GCMHD+ uses the same time varying A.C. as 
GCD+ (iKawata et al.ll2Qlll l, which reduces to zero where it 
is not required. Applying a minimum value of A.C. will pro- 
duce unnecessary smoothing in areas where it is not needed. 

Fig. I A2 1 shows the density profile produced by the code 
without it's scheme for A.C The wall heating error becomes 
significantly worse when the A.C. scheme is turned off. The 
wall heating problem shown in test 3A can only be solved 
by introducing a new switch for the A.C. which increases 
the strength of the dissipation when this type of shock front 
is detected. 

In Section 3.5 we showed how the V • B error varied 
as a function of the parameters OL'^in $ ^ind optimised 
the values of these parameters. The accuracy of the solution 
produced by the code changed when the parameters varied. 
This is easily seen in shocktube test 5A, as shown in Fig. 
1181 The V • B error reduces as the parameter /S is reduced. 
However, Figs. IA3llA5l show the effect of varying the param- 



eter $ on the solution produced for Vy. With = 0.2 the 
tensile instability is clearly not suppressed in the low reso- 
lution region, right of x = 0.15. The solution produced with 
/3 = 1.0 shows a deviation from the ATHENA solution at 
X = 0.15, where the velocity overshoots the reference solu- 
tion. In order to produce a Vx profile which agrees with the 
reference solution and has the least error a value of $ = 0.5 
was chosen. 



APPENDIX B: LATE TIME ORSZAG-TANG 
VORTEX 

In Section 3.3 the results for the Orszag-Tang vortex at t = 
0.5 were shown. At later times the test develops turbulence 
which then decays away. This development and subsequent 
decay of turbulence is very challenging for SPMHD codes 
to capture. We ran the Orszag-Tang vortex test again with 
a resolution of 600 x 692 and 1800 x 2076 and compared 
with the result produced by ATHENA with a resolution of 
600 X 600 cells. 

The result is shown in Fig. IBll The code correctly cap- 
tures the majority of the features present in the density field. 
However, the code fails to capture the density peak at the 
center of the simulation. As a result the central density fea- 
tures seen in the ATHENA simulations are not fully repro- 
duced in the GCMHD+ simulations. 



APPENDIX C: KELVIN-HELMHOLTZ 
INSTABILITY 

SPH codes in general struggl e to resolv e the Kelvin- 
Helmholtz instability (KHI). Ka wata et aD ([2009., ) demon- 
strated that the instability could be resolved when imple- 
mented with A.C We consider a periodic two dimensional 
box from x = —0.5 to 0.5 and from y = —0.5 to 0.5. Equal 
mass particles are set out on a square lattice with 724 par- 
ticles along the x-axis in the high density region, between 
y = —0.25 and y = 0.25, and 512 particles along the x-axis 
in the low density region. The high density region has a den- 
sity of ph = 1.0 and a velocity Vx = —0.5. The low density 
region has a density of pi — 0.5 and a velocity Vx — 0.5. The 
two regions are in pressure equilibrium with P/^ = = 2.5. 
The instability is seeded by adding random perturbations to 
the X and y components of the velocity with an amplitude 
of 0.01. 

Fig. ICll shows the result of the KHI test for the 
ATHENA reference code, GCMHD+ with a^^^ = 0.1 and 
GCMHD+ with a^^^ = 0.6. Two different GCMHD+ solu- 
tions are shown to show the effect of increasing the min- 
imum applied artificial viscosity. The instability develops 
more quickly in the grid method solution of ATHENA com- 
pared to the SPH solutions. This is clearly seen at both 
time steps displayed. The increase in minimum A.V. of 
GCMHD+ reduces the development of the instability con- 
siderably when compared to the lower GCMHD+ implemen- 
tation of Oi^n =0.1. 

We then performed an MHD KHI test. The same set up 
as above was used with the addition of a magnetic field in the 
X direction. A uniform and homogeneous magnetic field of 
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Figure Bl. Plot of the Orszag-Tang vortex at times t = 1.0 (upper row) and t = 2.5 (lower row) for ATHENA (left column), GCMHD+ 
[600 X 692] (middle column) and GCMHD+ [1800 x 2076] (right column). The arrows show the strength and direction of the magnetic 
field. The code fails to capture the central density feature even at the higher resolution and both resolutions are not perfectly symmetric. 




Figure CI. KHI test results for ATHENA (left), GCMHD+ [a^^n = ^-l] (middle) and GCMHD+ [a^^n = ^-^l (right) at t = 1.0 
(upper row) and t = 5.0 (lower row). The arrows show the velocity field. The superior ability of the grid code to resolve the KHI is clear 
when the ATHENA result is compared to the SPH results. 
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Figure C2. MHD KHI test results for ATHENA (left), GCMHD+ [a^^^ = 0.1] (middle) and GCMHD+ [n^^^ = 0.6] (right) at t = 1.0 
(upper row) and t = 5.0 (lower row). The presence of a magnetic field in the x direction prevents the growth of the instability in the 
GCMHD+ solutions. 

strength Bx = 0.129 was applied to all particles in the sim- 
ulation. Fig [C2l shows the results of this test for ATHENA, 
GCMHD+ (a^X^ = 0.1) and GCMHD+ (a^X^ = 0.6). 
The magnetic field stabilises the instability and prevents 
it from fully developing. The ATHENA result again shows 
greater development of the instability when compared to 
the SPMHD solutions, but the effect of the magnetic field is 
clear when the results are compared to Fig. ICll The solu- 
tions produced by the code show very minimal development 
of the instability. Further work is required to ensure that 
GCMHD+ is capable of resolving Kelvin-Helmholtz insta- 
bilities. 



